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ABSTRACT 

Very long baseline interferometry (VLBI) presents an opportunity to measure at the centimeter 
level such geodetic parameters as baseline length and instantaneous pole position. In order to 
achieve such precision, the geophysical and astronomical models used in data analysis must be as 
accurate as possible. The Mark-Ill interactive data analysis system includes a number of refinements 
beyond conventional practice in modeling precession, nutation, diurnal polar motion, UTl, solid 
earth tides, relativistic light deflection, and reduction to solar system barycentric coordinates. The 
algorithms and their effects on the recovered geodetic, geophysical, and astrometric parameters will 
be discussed. 
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INTRODUCTION 

If the earth were an isolated, rigid, homogeneous, spherical, airless body, the calculation of the 
VLBI observables, delay and delay rate, would be quite straightforward. In this simple system, the 
VLBI delay is given by the projection of the baseline, measured in units of light travel time, in the 
direction of the observed source; i.e., by the dot product of the baseline vector and the source unit 
vector. The delay rate is the derivative of delay with respect to time. On the idealized earth, the 
baseline is attached rigidly to the terrestrial coordinate system. Viewed from the celestial reference 
frame in which the sources are fixed, the components of the baseline vector change as the baseline 
is carried around by the rotation of the earth; consequently, the VLBI delay has a diurnal sinusoidal 
signature. Unfortunately for simplicity but perhaps fortunately for scientific interest, the situation 
of the real earth is more complex and a number of models are required to describe it. The object of 
most of these models is to compute the position and velocity of the stations and baseline in the 
celestial reference frame at the epoch of observation. This paper describes the models and algor- 
ithms which are used by the East Coast VLBI group and which are implicit in the work described by 
Robertson (1979) and by Knight (1979). 

It should be pointed out that other models for such phenomena as the propagation medium, an- 
tenna geometry, source structure, and clock instability are necessary for correct computation of 
the VLBI observables. As these effects will be considered in detail by other papers, they will not be 
discussed further here. 


Classes of Models 

As shown in figure 1, the models can be divided schematically into four classes. It should be em- 
phasized that this division is only schematic; the physics underlying different models is often closely 
related. The divisions follow the different modes of mathematical expression used. Models related 
to the rotation of the entire earth do not affect the baseline length but do affect its orientation and 
the orientation of the earth with respect to the celestial reference frame. Models related to the dis- 
placement of individual stations, such as solid earth tides and ocean loading, affect the baseline 
length. The model related to displacement of individual sources takes into account the gravitational 
deflection of the incoming radio signals by the sun. Models related to the orbital motion of the 
earth will not be discussed further since they are the subject of the paper to be given by Prof. 
Shapiro (1979). It should be noted, however, that the position and velocity of the earth with re- 
spect to the solar system barycenter are computed directly from the Planetary Ephemeris Program 
(PEP) tape provided by the Massachusetts Institute of Technology (MIT); consequently, no circu- 
lar approximations appear in aberration or potential. 

Figure 2 shows the mathematical form used for each class of model. The models in the first class 
are implemented as rotation matrices or products of rotation matrices; long period polar motion 
or wobble (W), diurnal polar motion (D), diurnal spin (S), nutation (N), and precession (P). These 
transform the components of the baseline from the terrestrial coordinate system to the reference 
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1. ROTATION OF ENTIRE EARTH 

2. DISPLACEMENT OF INDIVIDUAL STATION 

3. DISPLACEMENT OF INDIVIDUAL SOURCE 

4. ORBITAL MOTION 
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Figure 1 . Classes of models. 


Figure 2. Model algorithms. 


celestial coordinate system. The models of the second class are implemented as small translations 
of the station coordinates that vary with time for earth tides (ET) and ocean loading (OL). In fig- 
ure 2, the parentheses surrounding the second term on the second line indicate that this correction 
has not yet been used in our analyses. The single model of the third class is implemented as a small 
rotation of the source unit vector for gravitational deflection (GD). 


Rotation Models 

The purpose of the rotation models is to transform from the terrestrial coordinate system, in which 
the stations are initially located and in which, except for tectonic processes and short period tidal 
effects, the station coordinates are invariant, to the reference celestial coordinate system in which 
the VLBI observables are calculated. The origin of the terrestrial system is defined by the adopted 
coordinates of the intersection of antenna axes of the 37-m telescope at the Haystack Observatory 
(x = +1492406.691 m, y = -4457267.330 m, z = +4296882.102 m). These values were derived 
from the coordinates of the Mars antenna in the Goldstone complex of the Deep Space Network as 
determined from tracking data and the vector baseline from Goldstone to Haystack as determined 
from VLBI data. The reference celestial coordinate system is defined by the mean equator and 
equinox of 1950.0 with the origin at the solar system barycenter. The conceptual basis for the 
various rotations is that they represent the motion of the earth in space and the motion of the pole 
on the earth since 1950.0. 

Figure 3 shows the north polar region. The wobble matrix rotates the z-axis successively about the 
X- and y-axes from the geographic origin at 0 to the position of the slowly moving pole defined con- 
ventionally by optional instruments. Since our analyses have usually been performed considerably 
later than the VLBI observations were made, we use as a priori values the pole positions distributed 
by the Bureau International de I’Heure (BIH) in circular D. These initial pole positions are not 
necessary, and we normally recover polar motion when analyzing more than 1 day’s data. The 
second set of rotations for diurnal polar motion takes the z-axis from the slowly moving pole to the 
instantaneous pole. We use the seven largest terms in the series for the diurnal pole offset derived 
by McClure (1973), which has a maximum radius of 60 cm. While the effect is quite small (<.02 
arcsec), we feel that we have good evidence for it in our data. See C. Ma (1978) for a detailed 
discussion. 
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Figure 3. Wobble - long period polar motion 
diurnal polar motion. 

The matrix S for diurnal spin rotates the coordinate system about the z-axis in order to re-align the 
x-axis from the terrestrial origin of longitude at the Greenwich meridian to the celestial origin of 
right ascension at the vernal equinox. The rotation angle is the sum of four terms: (1) a constant 
sidereal rate from the classical equation of Newcomb given in the Explanatory Supplement to the 
American Ephemeris and Nautical Almanac (ESAENA), (2) the equation of the equinoxes for the 
difference between mean and apparent sidereal time, (3) an offset in UTl for variation in true rota- 
tion rate, and (4) four terms suggested by Woolard (1959) for theoretical short period variations in 
rotation rate. The a priori values for UTl offset are taken from circular D. As in the case of long 
period polar motion, the initial value for UTl offset is not necessary and is normally estimated if 
more than 1 day is analyzed. The terms suggested by Woolard are included because the smoothing 
used by the BIH removes short period variations. These terms have fortnightly and monthly periods 
with magnitudes less than 0.8 msec. The algorithm used to compute the Woolard terms was pro- 
vided by J. G. WUliams (1974).* 

The geometry for nutation is shown in figure 4. Instead of the commonly used first-order matrix 
for nutation, we apply an exact product of three rotations. The first rotation is about the x-axis 
(perpendicular to the flat figure at T) from the true equator to the mean ecliptic through an angle 
equal to the true obliquity ; i.e., the sum of the mean obliquity, e, and the nutation in obliquity, Ae. 
The second rotation is about the z-axis along the mean ecliptic by an angle equal to the nutation in 
longitude, A\]j. The last rotation is about the x-axis at M from the mean ecliptic to the mean equa- 
tor of date by an angle equal to the mean obliquity. The values for the nutations in obliquity and 
longitude come from the standard Woolard series given in the ESAENA as interpolated from tabu- 
lated values on the PEP ephemeris tape. In addition, seven terms suggested by Melchior (1971) to 
model elastic effects are used. The periods range from 18.6 years to fortnightly with amplitudes up 
to 0.02 arcsec in obliquity and 0.04 arcsec in longitude. The algorithm using the coefficients from 
Melchior’s table IXB was also provided by J. G. Williams.* 


♦J. G. Williams, JPL EM 391-592, private communication, 1974. 
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Figure 4. Nutation . 


It should be mentioned that the new nutation series to be proposed by the Working Group on Nuta- 
tion to the lAU in Montreal will affect the models for diurnal polar motion, diurnal rotation, and 
nutation. 

The model for precession follows conventional practice given in the ESAENA with rotations by z, 
0, and f. However, these angles are calculated from the constant of general precession suggested by 
Fricke (1967, 1971), which is 1.10 arcsec/century greater than the standard value. The algorithm 
allows the precession constant to be recovered, but the short span of VLBI data limits the precision 
of the estimate. The model and a priori source positions will require changes when epoch 2000.0 
coordinates are adopted. 


Station Displacement Models 

The tidal potential is calculated exactly using the equation shown in figure 5. The geocentric sta- 
tion vector is A, and the geocentric vector to the moon or sun, derived from the PEP ephemeris 
tape, is R. The distance from the station P to the attracting body M is p. 0 is the origin of the ter- 
testrial coordinate system. This total tidal potential includes both the second degree harmonic as 
well as higher degree harmonics, which are smaller by factors of A/R. To scale the station displace- 
ment from the tidal potential, we use the Love numbers calculated by Dahlen (1976; h = 0.609, 
1 = 0.085. Because of the geometry of the baselines we have used so far, the effect of solid earth 
tides on baseline length is only at the level of a few centimeters while the effect on baseline orienta- 
tion is negligible. The earth tide model has been used to estimate the Love numbers and a tidal lag 
angle. 

The aspect of ocean loading that affects VLBI observations is the displacement of a station as the 
tide moves masses of water from place to place, causing the crust to move in response. We have 
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Figure 5. Solid earth tides. 


investigated the model of Farrell (1977)* for the M2 and 01 tides. According to the model, the dis- 
placement at Haystack is less than 5 cm. Although we have prepared the computer code for imple- 
menting the model and scaling other tidal components, we have not applied ocean loading because 
of uncertainties as to its usefulness. 


Source Displacement Model 

The apparent displacement of a source caused by gravitational deflection of the incoming signal by 
the sun is modeled using a simple expression given by Shapiro (1967). The expression is applicable 
at all angles and allows the relativistic parameter gamma to be recovered. The necessary geometry 
between the source, the sun, and the earth is computed from the PEP ephemeris tape. 


Summary 

We have spent considerable effort in examining models, selecting and coding algorithms, and check- 
ing the results both internally and with other programs. Further details of these and other models 
can be found in Robertson (1977), C. Ma (1978), and the documentation of the CALC program in 
the Mark III data analysis system. Based on work done with older VLBI data, we believe that our 
geophysical and astronomical models are ready for the better data possible with the Mark III system. 

There are several areas in which we expect to improve our models, especially in ocean loading and 
nutation. We will examine other ocean loading models with different ocean tides and response func- 
tions. We will also investigate further the use of tidal gravimeter data from which the displacement 
might be scaled directly for the period of an experiment. For nutation, we are considering the for- 
mulation and series recently computed by Wahr (1979).** In addition, we are working with G. 
Kaplan of the U.S. Naval Observatory to estimate the coefficients of the nutation series using both 
VLBI and connected-link interferometer data. 

We expect that scientifically interesting results will continue to be forthcoming as the data allow the 
models to improve and as the models permit the data to be better understood. 


*W. E. Farrell, private communication, 1977. 

**J. N. Wahr, Ph.D. thesis, U. of Colorado, 1979. 
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